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Abstract 

We study cluster perturbation theory [Phys. Rev. Lett. 
84, 522 (2000)] when auxiliary field quantum Monte Carlo 
method is used for solving the cluster hamiltonian. As a 
case study, we calculate the spectral functions of the Hub- 
bard model in one and two dimensions and compare our 
results for the spectral functions to results obtained using 
exact diagonalization to solve the cluster hamiltonian. The 
main advantage of using quantum Monte Carlo results as a 
starting point is that the initial cluster size can be taken to 
be considerably larger and hence potentially capture more 
of the relevant physics. The drawback is that quantum 
Monte Carlo methods yield results at imaginary times with 
stochastic errors. 



1 Introduction 

In the field of strongly-correlated systems (such as 
the high-temperature superconductors), it is usually use- 
ful to study the single-particle excitation spectrum. Ex- 
perimentally, it can be obtained through the ARPES tech- 
nique for some systems |3|; while theoretically, it can 
be calculated by exact diagonalization (ED) or quantum 
Monte Carlo (QMC) method for some finite lattice mod- 
els, such as the Hubbard model, which describes a lattice 
system of strongly-interacting electrons and which is be- 
lieved by some Physicists to capture the Physics of the high- 
temperature superconductors. 

Since these calculations are for a finite system, it is desir- 
able to extend the calculations to the infinite lattice systems 
for a better comparison between the experiments and the 
theoretical calculations. Cluster perturbation theory (CPT) 
as proposed by Senechal et al. 1 1 01 1 1 II is one of these ex- 
tensions, and has turned out to be an accurate and econom- 
ical technique for calculating the spectral functions of the 
Hubbard model in one- and two-dimensional (ID and 2D) 
lattice systems. The method starts by dividing the infinite 
lattice system into a periodic array of clusters. The Hub- 
bard model inside a cluster is then solved by ED, and the 
hopping matrix between two neighboring clusters is treated 



as a perturbation. The cluster Green's functions from ED 
are then used to perturbatively construct the Green's func- 
tions of the infinite lattice. The method is economical in 
that one usually only diagonalizes a small cluster (typically 
less than 14 sites for the Hubbard model), and the infinite 
lattice Green's function with an arbitrary momentum k can 
then be constructed with a very small computational over- 
head. The computationally limiting step is therefore the ex- 
act diagonalization of the cluster which effectively limits 
its size to less than 20 sites or so for the Hubbard model. 
(The computer memory requirement grows exponentially 
with increasing cluster sizes.) On the other hand many ma- 
terials, such as molecular solids, display a "natural" cluster 
size (the molecule) which may be significantly beyond what 
can be treated with exact diagonalization techniques and it 
therefore becomes of great interest to study the accuracy of 
cluster perturbation theory when the cluster itself is treated 
with approximate techniques. 

In this paper we study the accuracy of the CPT 
method 1 1 01 1111 when auxiliary field quantum Monte Carlo 
(AFQMC) methods 0DI21II2I are used to solve the cluster 
Hamiltonian. To distinguish between these two methods, 
we call the usual CPT method using exact diagonalization 
methods to solve the cluster hamiltonian EDCPT and analo- 
gously the present method QMCPT. The advantage of QM- 
CPT is that we can deal with clusters as large as 60 sites (re- 
sults will be presented elsewhere) as long as the sign prob- 
lem 1 6 1, which is the appearance of negative probabilities in 
the QMC simulations, is not too severe. However, AFQMC 
yields imaginary time Green's functions with stochastic er- 
rors. In order to use such Green's functions as input to the 
zero temperature cluster perturbation theory formalism an 
analytic continuation to real times (frequency) has to be per- 
formed, a notoriously difficult step. Hence, it is not obvious 
that this approach will yield reliable results. Fortunately, as 
we show in the following, if high precision numerical data 
are available for the imaginary time Green's functions an 
analytical continuation using maximum entropy methods of 
the AFQMC data yield results that are in quite good agree- 
ment with EDCPT for the same cluster size. 

In the next section we briefly introduce the EDCPT 
method. This is followed by the introduction of AFQMC 



Figure 1. 1D and 2D clusters for the CPT cal- 
culations. In 1D a cluster of 12 sites is used 
in the paper. A 3 x 4 cluster is drawn for the 
2D square lattice. 
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Figure 2. Evolution of cluster Green's func- 
tion with imaginary time r in QMC simulation. 
We see that enough amount of data has made 
the error bars very small, which is essential 
in MEM |D. 



and QMCPT, describing details in the calculation. Spectral 
functions of the ID and 2D Hubbard Hamiltonians are then 
presented and compared to those obtained using EDCPT. 

2 Model 

The model Hamiltonian we consider is given by 



H = H + V, 



(1) 

(2) 



where 



(tj)a i 



is the Hubbard model inside the /'th cluster, and 



V 



( c n* c Jj° + h - c -) 



(4) 



(Ii,Jj)c 



is the hopping terms between two nearest neighbor (NN) 
sites i and j in two NN clusters / and J, respectively. Here 
the operator notations are standard for the usual Hubbard 
model, i.e., ^(cur) is the electron creation (annihilation) 
operator for spin a, and m a is the electron number opera- 
tor for site % with spin a. Inside the cluster U is the on-site 
Coulomb interaction energy, and t is the hopping integral 
between two NN sites. When t = t, we recover the stan- 
dard Hubbard model. In the following calculations t is used 
as the energy unit. 



3 EDCPT Formalism 

To compare EDCPT with QMCPT, we will briefly de- 
scribe the calculational steps of EDCPT for the Hubbard 
model. EDCPT first uses the Lanczos algorithm |5 1, which 
is an efficient method for obtaining several lowest eigen val- 
ues of sparse matrices, to diagonalize the Hubbard Hamil- 
tonian Hq of a cluster for its ground state energy Eq and the 
corresponding eigen state |^). The single-particle Green's 
function is then given by 



G°» = GL(w)+G&(w), 
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where e and h represent respectively the electron and hole 
part of the green's function, and 77 is a small positive param- 
eter to give the Lorenzian broadening of the delta functions. 
The calculation of these cluster Green's functions is stan- 
dard, and can be found in the literature (e.g., Ref. |9|). 

With these Green's functions G^{u>), the Green's func- 
tions for the superlattice G(K, u>) is constructed through the 
strong coupling perturbation | lOl ll II : 



Gy(K,w) = 



1 - V(K)G°(u 



(8) 



where K is superlattice's momentum vector, and V(K) is 
the Fourier transform of Eq. @. V(K) is given by 



R 



o,R,KR 

ij e > 



(9) 



2 



QMCPT, 1D, M=12 
OBC, <n>=1, p=10/t, U=4t 
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EDCPT, 1D, M=12 
OBC, <n>=1,T=0, U=4t 
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Figure 3. Single particle spectral functions of 
the 1D Hubbard model from (a) QMCPT and 
(b) EDCPT with open boundary conditions. 
Both methods produce roughly the same 
single-particle excitation energies (u> - fi)/t 
and energy gap value, although the spec- 
tral height is different at some energy values. 
See text for the reason. 



QMCPT, 1D, M=12 
PBC, <n>=1, p=10/t, U=4t 
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EDCPT, 1D, M=12 
PBC, <n>=1,T=0, U=4t 
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where R represents position of the cluster in the superlat- 
tice. A residual Fourier transform is performed to construct 
the lattice Green's functions in terms of the k vectors of the 
original lattice 

M 

G CPT (k,c) = - £ G«(k,a;)e- ik -< r '- r '> ) (10) 

where M is the number of lattice sites in one cluster. The 
spectral function of the infinite lattice is then given by 

A(k,w) = -- lim ImG CPT (k,w). (11) 



Figure 4. Comparison of single particle spec- 
tral functions of a 1 D Hubbard model from (a) 
QMCPT and (b) EDCPT with periodic bound- 
ary conditions. Both methods give roughly 
the same single-particle excitation energies 
at most energy values, but the spurious exci- 
tations inside the energy gap shows that QM- 
CPT with periodic boundary conditions is in- 
accurate. 



4 QMCPT Formalism 

In QMCPT the imaginary-time cluster Green's functions 
G°(k, t) are calculated using the auxiliary field QMC tech- 
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QMCPT, 2D, OBC 
<n>=1, p=10/t, U=8t 
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Figure 5. Single particle spectral functions 
of a 2D Hubbard model from (a) QMCPT and 
(b) EDCPT with open boundary conditions. 
The cluster is of dimension 3x4. Both 
methods predict the same single-particle ex- 
citation energies and energy gap, though at 
some energy values the spectral height is dif- 
ferent. See text for the reason. 



nique j|6ll71 ll2l . They are then analytically continued to ex- 
tract the real frequency spectral functions A°(k, uj) through 
the Maximum Entropy method (MEM) |8 1: 



G u (k,r) = / duj 



J A°(k, 



-0UJ 



(12) 



where (3 = 1/fcgT is the inverse temperature in the QMC 
simulation, and r is the imaginary time. The above inte- 
gral in Eq. d!2i is discretized with a small Aw and summed 
over a set of A°(k,uj) values. A smaller Auj produces 
a smoother image for A°(k,uj), but different Auj usually 
gives results in agreement with each other. We refer the in- 
terested readers to Ref. [8 1 for the details of MEM. The re- 
sulting A°(k, uj) can be used to construct the real frequency 
cluster Green's functions via |4] 



G°(k,cj) = J dui 



, A°(k,v') 
uj + irj — uj' 



(13) 



After a cluster Fourier transform from the k vectors to clus- 
ter i,j indices, and again following Eqs. ©, (fTJ and 
(lilt , we can evaluate the spectral functions A(k, uj) for the 
infinite lattice. 

Next we discuss in detail the calculation of space-time 
Green's functions with QMC simulations. Since we are 
considering only the cluster Hamiltonian, we can neglect 
Eq. (0} and write Eq. Q as 



W, 



where 



and 



(14) 



(15) 



(16) 



representing, respectively, the kinetic and potential energies 
of the cluster. 

The partition function is then given by 



= Tr\{e- A < T+W \ 



(17) 



where (3 — AtL, and L is the number of time slices in the 
imaginary time direction. After tracing out the fermionic 
operators, we get 

Z = 5^JJdet[l + Si(a)B£_i(a)-»Bi(a)] 
- ^detO({(r}, M ) T detO({<r},/i) i . (18) 
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The Bi matrices are defined as 



e -ArT e W(l) 



-u 



for z,j NN, 
otherwise, 



(19) 
(20) 
(21) 



where fi is chemical potential, tanh (A/2) = 
tanh(Ar?7/4), <Ji(l) = ±1 is the auxiliary Ising spin 
coupled with the electrons at lattice site i and time 
t = (I — l)Ar, and a — ±1 corresponds to f or j in Eq. 

In calculating the space-time Green's functions, we form 
a matrix of dimension ML x ML after a complete sweep 
over the M cluster sites and L time slices, the inverse of 
which gives the desired cluster Green's functions EE) 



V 

/ G°(l,l) 
G°(2,l) 



V G°(L,1) 
where, for /i > 1%, 
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G%{h,h) = (ci(h)c](l 2 )) 

= [B^B^-i ■ ■ ■ Bj 2 +i(l + 

Bi 2 ■ ■ ■ BiBl ■ ■ ■ Bi 1+ i) (22) 

After a Fourier transform of the above Green's functions, 
we have Gl(l, 1) = G°(k,r), where r = (Z - 1)At with 
(1 < / < L). They are the input for the MEM algorithm in 
Eq. (I12> . Note that in the above definition of finite tempera- 
ture Green's functions in Eq. 1221 we used the same symbol 
as that for zero temperature in Eq. (|5}, since there is no dif- 
ference in the application of the CPT method in these two 
cases. 

5 Results 

5.1 ID Case 

In this section we will present QMCPT results on the ID 
and 2D Hubbard models. See Fig^for the division of the 
ID and 2D lattices into clusters and the geometry of the 
respective clusters. We choose 12-site clusters in both ID 
and 2D cases, so that we can compare the results with those 



available in the literature. The intra-cluster hopping integral 
t is used as the energy unit, and the on-site Coulomb inter- 
action U = At or U = 8t. The inter-cluster hopping integral 
t will be set to t, too. In QMC the inverse temperature is 
/3 = 10/i, which will produce results close to the ground 
state. The Lorenzian broadening parameter is -q = O.lt, 
which is of the same order of the QMC temperature. 



In QMC simulations we discretize the imaginary time 
P = 10/t into 200 slices, i.e., At = 0.05/* to make 
the Trotter error smaller than the statistical ones. The 12- 
site Hubbard ring assumes either open boundary conditions 
(OBC) or periodic boundary conditions (PBC). One thou- 
sand complete sweeps over the space-time lattice are per- 
formed to warm up the system. For each boundary case, we 
collect 51,000 sets of space-time Green's functions, each 
100 sets of which are averaged as 1 bin (total 510 bins) 
for the subsequent MEM analysis. A typical G°(k, r) ver- 
sus r is shown in Fig. [2] For generating a smooth figure 
for A (In, to), the discretization of the real frequency is set 
at Auj = 0.025* in MEM. Note that a smaller Aw does 
not change the resultant figure except for making the figure 
smoother with longer computation times. Fig.|5]compares 
the QMCPT result with the EDCPT result for a ID Hub- 
bard chain constructed from the 12-site cluster with open 
boundary conditions. We find that QMCPT result agree 
very well with the EDCPT one. They both have a Hubbard 
gap opening at uj — \i = and the same positions of spectral 
peaks, although the peak heights are different in some posi- 
tions. Hence, a precise determination of the quasi-particle 
weights is probably not feasible with QMCPT. We expect 
this to happen, because MEM always broadens the sharp 
peaks or blurs the fine details of the exact results. In spite 
of this problem, QMCPT still produces spectral functions 
that satisfy the usual sum rules |4|. 



There has been a suggestion of using periodic bound- 
ary condition (PBC) in the clusters |2|, which, in ED, can 
greatly reduce the dimension of the Hilbert space through 
the translational symmetry. These additional hopping terms 
are then subtracted from the perturbation V during the con- 
struction of the superlattice Green's functions. Comparison 
of the EDCPT results with open and periodic boundary con- 
ditions are made in Ref. 1 1 1 1, where both results agree in the 
sense that a Hubbard gap opens at half filling for both cases. 
Here we repeat the same calculations in the same systems 
with both QMCPT and EDCPT. The results are shown in 
both Fig. 0] and Fig.gJ We find that QMCPT calculations 
with periodic boundary conditions produce spurious single 
particle peaks around u> — ^ = 0. Therefore, in what fol- 
lows we use open boundary conditions in the QMCPT cal- 
culations. 
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5.2 2D Case 

In the 2D Hubbard system, we choose a cluster of 
M x M y = 3x4 (see Fig.^ with open boundary conditions 
for QMC simulation. We have collected 71,000 sets of 
space-time Green's functions, which are averaged to yield 
710 bins for MEM analysis. We find that in this U = 8t case 
less sets of data makes the MEM analysis difficult 1 8 1. The 
QMCPT result together with the EDCPT result are shown 
in Fig. |5j where the k vector scan is along T — M — X — T 
(Please see the figure illustration of this momentum space 
scan in Ref. 1101 '). We see that the QMCPT produces re- 
sults close to those of EDCPT. (See also Ref. (10|-) The 
main peaks found in EDCPT are all preserved in the QM- 
CPT calculation, which again shows that we can safely use 
QMCPT for the lattice spectral function calculations (in that 
QMCPT predicts the correct single-particle excitation ener- 
gies and energy gaps though with different spectral heights 
from EDCPT at some energy values). This is again because 
of the broadening of sharp peaks and the blurring of fine 
details in the MEM analysis. 

6 Conclusions 

In summary, we have presented our first effort in apply- 
ing CPT method using AFQMC to solve the cluster hamil- 
tonian. At low temperatures our QMC and QMCPT results 
agree very well with those from EDCPT [TO| H H • Com- 
pared with EDCPT, QMCPT works at finite and very low 
temperatures, and the cluster sizes can be much larger. We 
expect QMCPT to be a useful tool in calculating the spec- 
tral functions in not only the ID and 2D lattices but also 
some molecular solids, e.g., fullerene materials, where each 
molecule naturally defines a cluster, and ED of Hubbard 
model is not possible for molecules with more than 20 sites. 
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